# -*- coding:utf-8 -*-
"""
Created on 2009-10-28

@author: Qiu wenfeng
"""

#def del2(f, *varargin):
#    '''
#    DEL2 Discrete Laplacian.ref matlab
#       L = DEL2(U), when U is a matrix, is a discrete approximation of
#       0.25*del^2 u = (d^2u/dx^2 + d^2u/dy^2)/4.  The matrix L is the same
#       size as U, with each element equal to the difference between an 
#       element of U and the average of its four neighbors.
#    
#       L = DEL2(U), when U is an N-D array, returns an approximation of
#       (del^2 u)/2/n, where n is ndims(u).
#    
#       L = DEL2(U,H), where H is a scalar, uses H as the spacing between
#       points in each direction (H=1 by default).
#    
#       L = DEL2(U,HX,HY), when U is 2-D, uses the spacing specified by HX
#       and HY. If HX is a scalar, it gives the spacing between points in
#       the x-direction. If HX is a vector, it must be of length SIZE(U,2)
#       and specifies the x-coordinates of the points.  Similarly, if HY
#       is a scalar, it gives the spacing between points in the
#       y-direction. If HY is a vector, it must be of length SIZE(U,1) and
#       specifies the y-coordinates of the points.
#    
#       L = DEL2(U,HX,HY,HZ,...), when U is N-D, uses the spacing given by
#       HX, HY, HZ, etc. 
#    
#       Class support for input U:
#          float: double, single
#        
#        ??? not full test!take care use
#    '''
#    msg = ''
#    nin = len(varargin)+1
#    loc = []
#    
#    # Flag vector case and column vector case.
#    ndim = f.ndim
#    vflag = 0 
#    cflag = 0
#    if f.ndim == 2:
#        if f.shape[1]==1:
#            ndim=1
#            vflag=1
#            cflag=0
#        elif f.shape[0]==1: # Treat row vector as a column vector 
#            ndim=1
#            vflag=1
#            cflag=1
#            f = f.transpose()
#    indx = f.shape
#    
#    # Default step sizes: hx = hy = hz = 1
#    if nin == 1:    # del2(f)
#        for k in range(f.ndim):
#            loc.append(numpy.arange(1,indx[k]+1, dtype=float))
#    elif nin==2:    # del2(f,h)
#        # Expand scalar step size
#        if len(v[0])==1:
#            for k in range(f.ndim):
#                h = v[0]
#                loc.append(h*numpy.arange(1,indx[k]+1, dtype=float))
#        # Check for vector case
#        elif vflag:
#            loc[0] = v[0]
#        else:
#            msg = 'Invalid inputs to DEL2.'
#    elif f.ndim == len(v): # del2(f,hx,hy,hz,...)
#        # Swap 1 and 2 since x is the second dimension and y is the first.
#        loc = v
#        if ndim>1:
#            tmp = loc[0]
#            loc[0] = loc[1]
#            loc[1] = tmp
#        
#        # replace any scalar step-size with corresponding position vector
#        for k in range(f.ndim):
#            if len(loc[k])==1:
#              loc[k] = loc[k]* numpy.arange(1,indx[k]+1, dtype=float) 
#    else:
#        msg = 'Invalid inputs to DEL2.'
#    
#    if msg<>'':
#        raise ValueError, 'Jolly:del2:InvaildInput%s'%(msg)
#    # Loop over each dimension. Permute so that the del2 is always taken along
#    # the columns.
#    if ndim == 1:
#        perm = numpy.array([[1, 2]])
#    else:
#        perm = numpy.concatenate((numpy.arange(2,ndim+1, dtype=float),1))
#    
#    v = numpy.zeros(f.shape, f.dtype)
#    
#    for k in range(ndim):
#        n, p = f.shape
#        x = loc[k][:]
#        h = numpy.diff(x, n=1, axis=0)
#        g = numpy.zeros(f.shape, f.dtype) # case of singleton dimension
#        
#        # Take centered second differences on interior points to compute g/2
#        if n>2:
#            g[1:n, :] = (numpy.concatenate((f[2:n+1,:],f[1:n,:])) / 
#                h.take(numpy.arange(2,n, dtype=float)).take(numpy.ones((p,1),float),axis=1) 
#                - numpy.concatenate((f[1:n,:],f[0:n-1,:])) /
#                h.take(numpy.arange(0,n-1, dtype=float)).take(numpy.ones((p,1),float),axis=1))/ \
#                (h.take(numpy.arange(1,n, dtype=float)).take(numpy.ones((p,1),float),axis=1) +
#                 h.take(numpy.arange(0,n-1, dtype=float)).take(numpy.ones((p,1),float),axis=1) )
#        
#        # Linearly extrapolate second differences from interior
#        if n>3:
#            g[0,:] = g[1,:]*(h[0]+h[1])/h[1] - g[2,:]*h[0]/h[2]
#            g[n-1,:] = -g[n-1,:]*(h[n-2])/h[n-3] + g[n-2,:]*(h[n-2]+h[n-3])/h[n-3]
#        elif n==3:
#            g[0,:] = g[1,:]
#            g[n-1, :] = g[1,:]
#        else:
#            g[0, :]=0
#            g[n-1, :]=0
#        
#        if ndim==1:
#            v=v+g
#        else:
#            v=v+ipermute(g, numpy.array([k-1:ndim+1, 0:k]))
#        
#        # Set up for next pass through the loop
#        f = permute(f, perm)
